import numpy as np
import matplotlib.pyplot as plt
def VV(x):
global L, V
if x>-L and x<L:
return 0
else:
return V
def f(r, x, E):
global h_bar, m
psi=r[0]; phi=r[1]
fpsi=phi
fphi=(-2*m/h_bar**2)*psi*(E-VV(x))
return np.array([fpsi, fphi], float)
def RK4(r, x, E):
global h
k1=h*f(r, x, E)
k2=h*f(r+0.5*k1, x+0.5*h, E)
k3=h*f(r+0.5*k2, x+0.5*h, E)
k4=h*f(r+k3, x+h, E)
r+=(k1+2*k2+2*k3+k4)/6
return r
def solve(E, alpha):
global L
psi=alpha; phi=1.
r=np.array([psi, phi], float)
for x in np.linspace(-1.5*L, 1.5*L, N):
r=RK4(r, x, E)
return r[0]
if __name__=="__main__":
L=1; V=20; h_bar=1; m=1
e=1.6022e-19
N=1000
x=np.linspace(-1.5*L, 1.5*L, N)
alpha=0.1
h=x[1]-x[0]
E1=0; E2=5
tolerence=1e-3
p1=solve(E1, alpha)
p2=solve(E2, alpha)
EE=(E1+E2)/2
while abs(solve(E2, alpha)-alpha)>tolerence:
if( (solve(E1, alpha)-alpha)*(solve(EE, alpha)-alpha)<=0):
E2=EE
else:
E1=EE
EE=(E1+E2)/2
ppsi1=np.zeros_like(x)
pphi1=np.zeros_like(x)
ppsi1[0]=alpha
pphi1[0]=1
for i in range(len(x)-1):
r=np.array([ppsi1[i], pphi1[i]], float)
r=RK4(r, x[i], E2)
ppsi1[i+1]=r[0]
pphi1[i+1]=r[1]
integ=0.
for i in range(len(x)):
integ+=h*ppsi1[i]**2
norm_ppsi1=ppsi1/np.sqrt(integ)
E1=3*EE; E2=10
p1=solve(E1, alpha)
p2=solve(E2, alpha)
EE=(E1+E2)/2
while abs(solve(E2, alpha)+alpha)>tolerence:
if( (solve(E1, alpha)+alpha)*(solve(EE, alpha)+alpha)<=0):
E2=EE
else:
E1=EE
EE=(E1+E2)/2
ppsi2=np.zeros_like(x)
pphi2=np.zeros_like(x)
ppsi2[0]=alpha
pphi2[0]=1
for i in range(len(x)-1):
r=np.array([ppsi2[i], pphi2[i]], float)
r=RK4(r, x[i], E2)
ppsi2[i+1]=r[0]
pphi2[i+1]=r[1]
integ=0.
for i in range(len(x)):
integ+=h*abs(ppsi2[i])**2
norm_ppsi2=ppsi2/np.sqrt(integ)
plt.plot(x, norm_ppsi1)
plt.plot(x, norm_ppsi2)
plt.show()